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Abstract. We consider the A-body problem in a layered geometry containing cold polar molecules with 
dipole moments that are polarized perpendicular to the layers. A harmonic approximation is used to 
simplify the Hamiltonian and bound state properties of the two-body inter-layer dipolar potential are used 
to adjust this effective interaction. To model the intra-layer repulsion of the polar molecules, we introduce a 
repulsive inter-molecule harmonic potential and discuss how its strength can be related to the real dipolar 
potential. However, to explore different structures with more than one molecule in each layer, we treat 
the repulsive harmonic strength as an independent variable in the problem. Single chains containing one 
molecule in each layer, as well as multi-chain structures in many layers are discussed and their energies 
and radii determined. We extract the normal modes of the various systems as measures of their volatility 
and eventually of instability, and compare our findings to the excitations in crystals. We find modes that 
can be classified as either chains vibrating in phase or as layers vibrating against each other. The former 
correspond to acoustic and the latter to optical phonons. For the acoustic modes, our model predicts a 
smaller sound speed than one would naively get from expansion of the dipolar potential to second order 
around the origin. Instabilities can occur for large intra-layer repulsion and produce diverging amplitudes 
of molecules in the outer layers, and our model predicts how the breakup takes places. Lastly, we consider 
experimentally relevant regimes to observe the structures. The harmonic model considerd here predicts 
that for the multi-layer systems under current study chains with one molecule in each layer are always 
bound whereas two chains comprised of two molecules in each layer will not be bound. However, since 
realistic systems have external confinement prevention the molecules from escaping to infinity, we still 
expect the unstable modes to show up as resonances in the dynamics. 

PACS. 67.85.-d Ultracold gases, trapped gases - 36.20.-r Macromolecules and polymer molecules - 
03.75.Kk Dynamic properties of condensates 



1 Introduction 

The experimental study of cold dipolar molecules is a 
rapidly accelerating field [HHISlllISlilllllSlinilin] ■ The long- 
range anisotropic forces of polar molecules can be con- 
trolled by external alignment and tuned to be both re- 
pulsive and attractive depending on the geometric setup, 
which has lead to a significant amount of interesting theo- 
retical proposals pTlfT^ . In particular, stacks of thin lay- 
ers containing polarized molecules interacting via dipole- 
dipole potentials have been suggested |!T3 l fT4 l [T5 p i6 | IT7] as 
a way to control the losses that can severely influence 
three-dimensional experiments with strong dipolar forces 
[TllStlTSlfTQ] . The layered structure and the long-range na- 
ture of the interactions holds promise for the realization of 
interesting few- [20,21,22,23,24,26,27,28,29,30] and many- 
body states [3lll^2t.33,,34t.35^.36,.37,,,38,.39t.40..41t.42). Recently,?^!^"^^*!™^ [60,61^. 



tallization when the dipole moment becomes large. This 
phenomena is similar to the famous Wigner crystal phase 
of the electron gas [43ll44]. However, whereas the electron 
gas crystallizes at low density where the Coulomb inter- 
action dominates, a system of dipoles will crystallize at 
high density where the kinetic energy is negligible com- 
pared to the dipole force. In a single two-dimensional (2D) 
layer with dipoles oriented perpendicular to the layer this 
has been studied in classical simulations , and more 

recently in quantum Monte Carlo simulations [THHSIHS] . 
and for large dipole moments evidence for a triangular 
crystal structure is found. Studies in one-dimensional (ID) 
dipolar systems also find interesting crystal phases [SCTIISTI 
[52 | f53 l [54 l l55 p 56| , l57] which are similar to those seen in ion 
Coulomb crystals [551155] . Crystal phases have also been 
found in 2D with an in-plane optical lattice for arbitrary 



a stack of layers have been realized with fermionic polar 
^•^K^^Rb molecules [1^ and further cooling should pro- 
duce some of the novel phases that have been proposed. 

An interesting question to pose for layered systems of 
dipolar molecules is related to their tendency for crys- 



Here we are interested in the case of a multi-layer sys- 
tem with dipoles perpendicular to the planes as in re- 
cent experiments [lOj . The possibility of having bound 
states consisting of chains with one molecule in each layer 
has been discussed for both bosonic [T7ll62l and fermionic 
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molecules |63j . In the limit of large dipole moments, a 
study of the two-layer (bilayer) case using classical dy- 
namics showed that the triangular crystals appear spa- 
tially correlated, i.e. the molecules in the crystal lattice 
sit on top of each other [51] . This is also the expected be- 
haviour for more than two layers and we note that a recent 
quantum Monte Carlo study finds evidence of a transition 
to rough chains for bosonic molecules at a critical dipole 
strength [65] , 

In the present work we consider strongly interacting 
distinguishable dipolar molecules before they enter the 
crystal phase. Rather than fixing the position of the molecule: 
to a particular lattice and calculate fluctuations about this 
state, we consider an A^-body system with a fixed number 
of molecules in each layer and study the quantum spec- 
trum to obtain information about the modes of the system 
before it crystallizes. Most of the studies on multi-layer 
systems mentioned above are based in one way or another 
on a harmonic approximation to the interaction potential 
which renders the otherwise intractable A'^-body problem 
solvable in certain limits, and provides ground state and 
excitation modes of the system. In this paper we also make 
use of the harmonic approximation to derive a solvable N- 
body problem. However, we go beyond previous studies in 
two important aspects; i) we use the two-body binding 
energy to fix the interlayer attractive interactions for the 
effective harmonic hamiltonian and ii) we include para- 
metrically the repulsive interaction between molecules in 
the same layer which has been neglected thus far except 
in the classical calculations of Ref . [53] . 

To arrive at our effective harmonic Hamiltonian, we 
use the method recently formulated in Refs. [271155] . We 
proceed in two steps, that is first the one- and two-body in- 
teractions are approximated by quadratic forms, and sec- 
ond the resulting A^-body Schrodinger equation is solved 
exactly. We replace the true interactions with quadratic 
forms in the molecule coordinates, either by direct fits of 
the potentials or by adjusting parameters to reproduce 
crucial properties. This approximation allows analytical 
investigations of the A"-body system with the properties 
expressed in terms of the two-body characteristics. 

The harmonic model studied here predicts that chains 
with a single molecule are stable for any interaction strength. 
A complex with two molecules in each layer, i.e. two chains 
in close proximity, will most likely not be bound as we es- 
timate the critical strength for breakup to be below the 
strength in current experimental setups. However, experi- 
ments are always performed in the presence of an external 
confinement and we expect the unstable modes of com- 
plexes with multiple chains to appear as resonances. We 
therefore study the modes as function of the intralayer 
repulsion to determine how the system breaks into single 
chains and elucidate its structure. 

The paper is organized as follows. In Section [2] we 
describe the harmonic approximation method used and 
Section [3] discusses in detail the chain structure with one 
molecule in each layer. This involves energies, radii, and 
normal mode excitations. Section [4] deals with multiple 
chains and layers with many molecules. The intralayer re- 




Fig. 1. Pictorial view of the multi-layered system containing 
dipolar molecules. The attractive interlayer (Va) and repulsive 
intralayer (Vr) are illustrated, along with a potential bound 
chain state containing four molecules. 



pulsion is taken as a parameter in the harmonic approx- 
imation. The normal modes reveal the most likely decay 
mechanism, and extract configurations of the most unsta- 
ble degrees of freedom. In particular, we calculate the crit- 
ical stability properties of the system and discuss how the 
system breaks up into smaller structures. In Section [5] we 
discuss how the intralayer repulsion influences the system 
size and densities, and we discuss relative energies of the 
many-body systems. Finally, Section [6] contains summary 
and conclusions. 



2 Method 

We consider a system of A^ dipolar molecules of mass, m, 
and dipole moment, D, distributed in a series of paral- 
lel two-dimensional layers separated by a distance d as 
illustrated in figure [T] The planar coordinates of the k'th 
molecule are {xk,yk)- The dipoles are oriented perpen- 
dicular to the layers and the molecules interact pairwise 
through the dipole-dipole potential, V , which in different 
planes is given as 



V{x,y,n)^D^ 



2{ndf 



(a;2-h2/2 + („(f)2)5/2 



(1) 



where n is the number of layers separating the two molecules, 
and (x, y) are the relative coordinates of the two molecules. 
We define the dimensionless strength, U — mD^ /{h^d). 
For n = 1 we have a nearest neighbour interaction, for n — 
2 a next-nearest neighbour interaction and so forth. Below 
we will also use the standard notation for the radius r = 
+ 2/2. Molecules in the same layer repel each other 
and if identical are also subject to (anti)symmetrization 
conditions for fermions and bosons. Here we will ignore 
symmetry considerations and consider our molecules dis- 
tinguishable. Since we work in the strongly-coupled regime 
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r/d 

Fig. 2. Potential energy curves for the nearest neighbour 
dipole-dipole interaction {Vdip{n — 1)) and the harmonic oscil- 
lator potential (Vhoin = 1)) with the same binding energy. The 
next-nearest neighbour dipole-dipole interaction is also shown 
{Vdip{n — 2)) as well as the oscillator potential with the match- 
ing binding energy {Vho(n ~ 2)). The dipole strength is here 
chosen as U — 6. The resulting binding energies are B2 ~ 3.44 
and 0.20 in units of /mtf . 




Fig. 3. Binding energies for two, B2, three B3, and four B4, 
molecules as functions of strength, U = niD^ /h^d, of the dipole 
potential. The oscillator approximation (solid, red) and the 
accurate numerical calculation (dashed, blue) are compared for 
three (_B3(num)) and four molecules (_B4(num)). The numerical 
calculations are courtesy of A. G. Volosniev using the method 
described in [2811291. 



and we do not expect the chains to have large overlap this 
is a reasonable assumption. The Hamiltonian for the N- 
body system is then 



H = -- 



N 

2m ^ 

k=l 



dxl 



dyl 



(2) 

where {xik,yik) = {xt - Xk,Vi - yk), and V for different 
layers is from equation ([l]) or the corresponding repulsion 
for molecules in the same layer. This A'^-body Hamilto- 
nian is solved using the method described in Ref. [66] . 
Note that there is no external trapping potential in any 
of the planes. The A^-body structures we consider are self- 
bound, i.e. they are held together by the attractive inter- 
actions between the layers. For configurations that have 
more than one molecules in a single layer, the attractive 
forces will have to overcome the repulsive intra-layer in- 
teractions. This will be disucssed in great detail later and 
will give raise to critical values of this intra-layer repulsion 
above which there are no bound structures. 

To render properties of systems with A^ molecules an- 
alytically tractable, we simulate the effects of the dipole- 
dipole interaction by harmonic oscillator potentials as de- 
scribed in "BB . We replace the actual interaction in equa- 
tion ([1]) with a shifted harmonic oscillator, Vtoi which has 
its node in the same place as the dipole-dipole potential. 



I.e., 



2{ndY 



(3) 



The potentials in equations ([T]) and (|3| are negative (and 
positive) in the same regions, and have nodes in the same 
places, see figure [2] The strength, Vbi is finally adjusted 



to reproduce the two-body dipole-dipole binding energy. 
This means that Vq is a specific function of the dipole mo- 
ment, the molecule mass, and distance between the layers, 
Vb(-D^,m, d) The corresponding oscillator frequency, Wa, 
is then given by 



2Vo{D'^,m,d) 



(4) 



where the reduced mass of the two-body system is m/2. 

The key ingredient is here the two-body binding en- 
ergy, B2, which from equation ^ for two molecules is 
seen to scale with layer distance, nd, as 



B (n) - f(^^^ 
m{nd)^ \ndlT? 



(5) 



where / is a universal function. For n = 1, that is, for 
neighbouring layers, we show B2 in figure [3] as function of 
potential strength, . The two-body interaction for pairs 
of molecules in layers separated by the distance, nd, is then 
obtained through the scaling in equation ([s]). Upon com- 
paring the root-mean-square radius of the exact solution 
to the oscillator approximation one finds very good agree- 
ment for U ^ 2 which is also where the binding energy of 
the exact and harmonic approximation for two molecules 
are very close [MIED]- 

These energies are then the crucial input in the A^- 
body calculations. At small strength the approach to zero 
is extremely strong [281129] whereas at larger strengths the 
behaviour is smooth and in a region almost linear in U. 
We need two-body interactions between each pair, the 
most important being between adjacent layers at a dis- 
tance d. For larger than one-layer separations the poten- 
tial becomes increasingly more shallow, flat and spatially 
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extended, see figure[2j The scaling of B2 with layer number 
is computed for a strength U/n (see equation ([5|). 

This procedure of adjusting strengths for different shapes 
of potentials to get the same binding energy has proved 
efficient in other contexts where effective interactions are 
used. The method is frequently applied in dripline nuclear 
physics, where knowledge about the systems can be lim- 
ited to an energy of one bound state, or resonance [7T]- It 
is also the philosophy employed for cold atoms where the 
scattering length is used as the parameter characterizing a 
short-range two-body interaction, and subsequently used 
in A'^-body calculations. However, the dipolar interaction 
is long-ranged. In fact, for typical experimental densities, 
the range of the dipole-dipole forces is comparable to the 
inter-particle distance. This implies that one should take 
more than just the zero-energy properties into account. 

The accuracy of the procedure is tested by comparing 
the oscillator approximation with numerical calculations 
of the energies of three and four molecules with one in 
each layer using a method similar to the one described 
in p51[2^ , see figure [sj The remarkable precision of the 
oscillator results for three molecules is reassuring in the 
applications on N molecules. 

We note that the procedure used here to fix the fre- 
quency of the attractive interlayer interaction is naturally 
very different from expanding the potential in equation 
([I]) for r/d <C 1 to get the harmonic term [30J. We include 
more accurately the two-body properties into the A^-body 
system and expect therefore to provide a better approxi- 
mation for the dynamics of a single chain and also for the 
interaction of several chains as we discuss below. 



2.1 Repulsive In-Plane Interactions 

Two polar molecules in the same layer repel each other 
when both dipoles are perpendicular to the plane. This re- 
pulsive intralayer interaction scales with the dipole strength 
similarly to the attractive interlayer dipolar interaction. 
Within the framework of effective harmonic hamiltonian 
we would like to replace the real intralayer potential by 
something that can model a repulsion, which can be achieved 
by using an inverted oscillator. The repulsive interaction 
is therefore replaced by a term of the form 



We will demonstrate below that the inverted oscillator 
provides us with a parametrically sound way of including 
repulsive forces in our model. That is, it provides a handle 
on the effect of repulsion in the detailed balance between 
repulsively and attractively interacting pairs of molecules. 
Therefore we are able to consider various properties of 
several chains in multiple layers which is one of the goals 
of this work. 

At this point we must caution that one might naively 
try to fit the frequency of the repulsive term above to 
the true intralayer dipole interaction which has the form 
V{r) = X/r^. Here we introduce the dipolar strength A 
since we want to separate it from the strength of the at- 
tractive interlayer potential (IT]) for later discussion. If we 
have a system at a certain density, n, then the charac- 
teristic energy of the true repulsion can be estimated as 
energy scale of an oscillator is usually hwr and 
we then arrive at 



mX 
Wd 



(8) 



1 2 

4 



(6) 



where 1/4 arise from use of the reduced mass, ujr is the cor- 
responding frequency, and Vro is a constant shift. This in- 
verted oscillator favors a large (infinite) distance between 
the molecules, and acts consequently as a repulsion. The 
shift is parametrized as 



-ma oj: 



2' -'d^ 



(7) 



where we use d as the natural length, and leave the dimen- 
sionless scale factor, a, for later adjustment. The distance 
where the potential changes sign is then r = ad. 



For current experiments |T0; with mX/h^d ~ 0.1 and nd^ ~ 
1, this gives a very small contribution from the repulsive 
term. However, the procedure of harmonic approximation 
breaks down for small U and/or A, and both the attractive 
and repulsive parts of the potential cannot be modelled by 
an oscillator in this limit. Only in the limit of large U does 
the approximation become valid for the single-chain struc- 
tures as shown above. Quantitatively, we expect that har- 
monic oscillator wave functions are reliable for [/ > 2 [571 
l30] . Below this value, the binding energy decreases rapidly 
and the state becomes extended [^[^ . 

As will be shown below, there are critical values of 
uJr{U) above which the system becomes unstable as indi- 
cated by the normal mode frequencies that become imag- 
inary. If we take densities nd^ ^ 1, the oj^ given by equa- 
tion Q is always larger than the critical frequency and 
the system is unstable. Furthermore, the densities we cal- 
culate below for the system are much larger than nd^ — 1, 
so this gives an apparent inconsistency. 

There are, however, two major points that invalidate 
the identification in equation ([8). First, there is no zero- 
point energy scale for an inverted oscillator, and using hujr 
as such is meaningless. One might attempt to make sure 
that the force (gradient of the potential) of the true poten- 
tial and the inverted oscillator match in a region around 
zero, but this suffers from the same inconsistencies. The 
second, and more important objection, is that we must 
choose the repulsion in a manner that is consistent with 
the philosophy of using effective harmonic approximations 
already employed for the attractive interlayer interaction 
discussed above. There we really use the knowledge of two- 
body states in the system to fix a potential that repro- 
duces features such as energy and potential shape. 

In order to apply our philosophy, we can compare to 
small systems that have the effect of the repulsion in- 
cluded, but which can be solved by other means. The 
obvious choice is the four-body state with two molecules 
in each of two layers or two molecules in one layer and 
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one molecule in each of the two adjacent layers in a setup 
with three layers. We could then compare the energetics 
and fit the strength of the repulsive term. Exact results 
on these configurations have been reported very recently 
in Ref. |25j. These results indicate that such systems are 
actually unstable for a critical intralayer repulsion that is 
much smaller than the situation where attraction and re- 
pulsion have the same value {U = X below). Ref. [25 also 
find that longer chains (five layers or more) could bind ex- 
tra molecules as long as the additional ones are not in ad- 
jacent layers. This hints at a competition between attrac- 
tive and repulsive interaction terms that is conveniently 
described within the harmonic oscillator approach as we 
demonstrate below. The results of Ref. [2S] are, however, 
difficult to use as a fit for the repulsive frequency in the 
harmonic model due to the instability of the states. We 
therefore follow a different line of argument which relates 
the number of bound states in the (inverted) repulsive 
potentials as we now explain. 

In a single plane, the pure inverse cubic dipolar repul- 
sion obviously has no two-body bound states. However, 
if we introduce an external trapping potential, i.e. an ex- 
ternal oscillator to keep the molecules confined, then the 
two-body spectrum would be shifted upward due to the 
repulsion. One could then imagine adding a repulsive oscil- 
lator term to fix the two parameters of the potential in (|6| 
so that, say, the two lowest bound state energies match (or 
the ground state and the root- mean-square radius) . Unfor- 
tunately, the results of such a procedure depend strongly 
on the choice of external oscillator potential. If one relates 
this external confinement to the density of particles in the 
layer (in a manner similar to ([s])), then we are left with 
a two-body interaction that depends on the properties of 
the entire system. We consider this situation extremely in- 
convenient. This dependence of two-body physics on the 
many-body problem is also inherent in the naive estimate 
in (g). 

What can be done instead, is to consider the num- 
ber of bound states that the inverse dipolar potential, i.e. 
an attractive potential, allows and compare this to —Vrep 
(|6| (a normal oscillator that is shifted below zero at the 
origin). We consider only those bound states that have 
negative energy even though the shifted normal oscillator 
will of course have the usual ladder spectrum of states. 
This explains the upper limit in (11 1 below. We are inter- 



ested in the stability of structures with several molecules 
in each layer, and therefore we need only require that the 
repulsion be reproduced within the range where there is 
an attraction to particles in other layers. This means that 
the shift term in ([t]) should ensure that the inverted oscil- 
lator crosses zero along the outer repulsive barrier of 0. 
This is roughly in the region a/2 < a < 2\/2 (the upper 
bound being defined as the point at which the repulsion 
of ([!]) is only 1% of its attractive value at the origin). The 
parameter a is therefore still somewhat arbitrary although 
we will compute a suitable value based on the energetics 
of barely stable configurations of molecules in section [5] 

In a realistic setup, the layers have finite width, w, and 
therefore the 1 /r^ repulsion is modified at small distances 



[26j and attains a logarithmic dependence on relative dis- 
tance of two molecules. The crossover happens for dis- 
tances of order the layer width (typically w ~ O.ld — 0.2d 
in experiments). Matching at w, we have 



Vdtp{r) 



-3Aln(- 



„l/3 , 



r > w 
r < w 



(9) 



An estimate for the number of bound states in a two- 
dimensional potential can be found by integrating the (ab- 
solute value of the) negative part of the potential over 
space [67^ . For the dipolar potential we find 



Vdipir)rdr 



9 A 

4w' 



while for the oscillator of ([6| we find 

! 1 
Vrep(''')r dr — — raiJld'^a^ 
16 



(10) 



(11) 



If we equate these expressions we obtain a relation be- 
tween the repulsive frequency, oj^, and the parameters A, 
a, and w, which reads 



6 



d I mX 
w V ti^d' 



(12) 



in dimensionless form. What we have done here is to use 
the number of bound states of the inverted potentials, 
which implies that we have in fact matched the number 
of bound states excluded by the presence of the repulsion. 
This is in similar spirit to the proposal of introducing an 
external oscillator discussed above. However, we avoid a 
dependence on properties beyond two-body physics. 

In section [4] we calculate critical values of tUr beyond 
which configuration with several molecules in each layer 
become unstable and break into smaller complexes. The 
critical repulsive frequencies we find below are md'^uir/h ~ 
1 — 10. Taking the specific case where mX/h^d = 20, we 
find md'^uJr/h = 7.78. From Q with w = 0.2d and 

= 2, we find a larger value of md'^ujr/h = 30. This 
implies that we are in the unstable regime. However, we 
note that the parameter a is geometrical and somewhat 
arbitrary, constrained only by our desire to describe the 
competition of repulsive and attractive terms in the sys- 
tem. Our calculations below indicate that a better value 
is a ~ 1.92. This yields md'^uJr/h ~ 16, which is closer to 
the critical frequency. 

We are thus in a situation where structures contain- 
ing more than one molecule in a single layer are most 
likely unstable in realistic experiments. This is consistent 
with the conclusion of Ref. [IS] about the instability small 
complexes. However, as noted above the experiments have 
external confinement and the structures obtained for mul- 
tiple chains below the critical frequency in later section 
could therefore appear as resonances in the confined sys- 
tem. The energetics and the path to breakup for these sys- 
tems is therefore an interesting question nonetheless. With 
a harmonic model we can study such questions for large 
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number of particles in an essentially exact manner, given 
the necessary approximation on the interaction terms. 

In making a distinction between mX/fPd and U (from 
the interlayer potential Q), we also leave open the pos- 
sibility of working with more than one type of molecule 
or populating different internal rotational states of the 
molecules in different layers, both of which can modify 
the dipolar moment and potentially bind the structures we 
study here by reducing the overall repulsion. This could 
possibly be achieved by externally applied electromagnetic 
fields and lasers [68,69,70 . However, in the present study 
we will always assume that X ^ U and use just U as 
the dipolar strength parameter. In light of the extended 
discussion above, we apply repulsive intralayer interaction 
parametrically through the inverted oscillator to study the 
effects of inter-chain interaction for the purpose of the 
present work. 

The oscillator potential is now defined for two molecules 
either in different layers or the same layer. If a confinement 
is needed it is straightforward to add and include a one- 
body harmonic oscillator potential. Without a one-body 
external potential the relative motion separates from the 
free motion of the center of mass which therefore becomes 
uninteresting. We shall therefore consider only the relative 
motion. 



3 One molecule per layer 

The simplest many-body structure is found when only one 
molecule is placed in each layer. This means that all the 
two-body interactions are attractive and given in equation 
([T]|, and the repulsion in equation (|6| is not present. Any 
pair of these molecules would then form a bound state 
as shown previously [66,28,29 . When the layers are far 
apart the potential is very shallow and the attractive part 
extends to a distance, r = ndy/2, proportional to the dis- 
tance between the layers. The attraction decreases with 
the third power of the layer distance, and the binding en- 
ergy decreases correspondingly Therefore the interactions 
dominate between the nearest neighbours. 



1.12 
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N 

Fig. 4. Energies divided by two-body energy obtained with 
the potential in equation ([T]) times number of pairs of nearest 
neighbours for several different dipole strengths. These energies 
were obtained by restricting the potential to nearest neighbours 
only. 
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Fig. 5. Energies divided by two-body energy obtained with 
the potential in equation ([T]) times number of pairs of nearest 
neighbours for several different dipole strengths. 



3.1 Binding energies and radii 

The energies could be expected to be proportional to the 
two-body energy, and linear in the number of inter- 
acting neighbours pairs, — 1. This was conjectured and 
estimated for large N in i.e. 



Bn 
{N-1)B2 



00 _^ 



TT 



1.645 



(13) 



However, the assumptions are then that the interactions 
decrease as l/n^, and that no additional correlations arise 
from the string of molecules. A better estimate should 
then be 



B 



N 



{N-1)B2 



00 - 

Y- 



C(3) w 1.202 



(14) 



where the first two powers are from the d^-dependence and 
one power is from the reduction of interaction strength, see 
equation ([5]). The latter is an upper limit since the energy 
curve in figure [3] is concave. With only nearest neighbour 
interactions we get the energies shown in figure [4] for dif- 
ferent dipole strengths. The variation is largest for small 
N while the energy per interacting pair is leveling out 
for large A^. All curves are between 1 and 1.12 with a 
monotonous decrease with interaction strength. Thus cor- 
relations always contribute on less than the 10% level, and 
increase as the system becomes spatially more extended 
and dilute. 

Inclusion of all interactions give the energies shown in 
figure [5] The trends are the same with increase as func- 
tion of A^ and decrease as function of interaction strength. 
However, the actual values are now all between 1 and 1.6. 
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Fig. 6. Energies in Fig [4] divided by the energies in figure 
[5] which is the fraction of the total energy of the A^-molecule 
string captured by the interaction of nearest neighbours as a 
function of chain length for several different dipole strengths. 
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Fig. 7. Radius-squared as a function of the position in the 
chain of 30 molecules. Position one is the outermost position 
of the chain, and the radii are symmetric with respect to the 
center of the chain, so showing only the first half of the chain 
is necessary. 



The large- TV asymptotics seems all to be above 1.202 in- 
dicating additional contributions from correlations which 
seems to approach zero in the limit of strong interactions. 

The 20-50% difference between energies in Figs, [l] and 
[5] is due to nearest neighbour and inclusion of all contribu- 
tions to the energy. The corresponding fraction is directly 
shown in figure [6] where the trend is decrease of nearest 
neighbour contribution with N and increase with interac- 
tion strength. Again the strongest variation is for small N 
while asymptotic constant ratios are approached for large 
A''. All ratios are between 0.7 and 0.9, showing how good 
an approximation only nearest neighbours would be, de- 
pending on strength and molecule number. The total num- 
ber of interaction terms, (TV — 1)^/2, increase compared to 
the number, A'^ — 1, of nearest neighbour terms. However, 
the rather strong fall off of the interaction with distance 
leads to the relatively small values in figure [6) Even when 
a chain gets longer, only the very nearest links in the chain 
are affected, and constant ratios are approached. 

The approximately linear dependence of energy on A^ 
implies that the separation energy, B]^_i — Bpf oc -B2, 
where the proportionality constant crudely can be approx- 
imated by C(3) = 1.202 (equation (14|) for large A^ and 



not too weak interactions. The variation with A^ of the 
separation energies is around 0.2 for all strengths. The 
correlations therefore change smoothly with the number 
of molecules. 

The wave functions for each of the states are also avail- 
able. They are essentially products of Gaussians of linear 
combinations of molecule coordinates. To get an idea of 
the structure, we computed mean square radii of the indi- 
vidual molecules measured from the common (projected) 
center of mass of all molecules. These radii are shown in 
figure [7] for a system of 30 molecules distributed in 30 lay- 
ers as function of the number of layers. The molecules in 
the outer layers have largest radii intuitively correlated 
with the fewer efficiently interacting neighbours. Towards 




Fig. 8. Radius-squared of the central molecule in a chain of 
length A as a function of chain length for several different 
dipole strengths. 



the central layer, the radii decrease and become almost 
independent of layer position. The radii decrease with the 
strength of the interaction in agreement with the smaller 
binding energies. This is consistent with the variational 
study in Ref. P71. 

The absolute sizes are on average between 0.3d and 
0.9d for the strengths in figure [7] Clearly d is the natural 
unit since that is the extension of the potential within each 
of the layers. When a radius is smaller than d it means that 
the molecule mostly is located inside the attractive well 
from the nearest neighbour. 

To see how the radii evolve with molecule number, we 
show the central radius as a function of molecule number 
in figure |8] For the same interaction strength, the size 
increases with molecule number. The outermost molecules 
are always the largest, as they only have the influence of 
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Fig. 9. Normal mode frequencies (points) as a function of 
molecule number for a dipole strength U — 6. The line is the 
average value of the frequency at a given number of molecules. 
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Fig. 10. Normal mode frequencies (points) as a function of 
number of molecules along with their fits to an A/N + B 
curve, where A and B are fitting parameters for a dipole 
strength (7 = 6. 



one nearest neighbour. The size of the outermost molecule 
also grows slightly faster than the central molecules. 



3.2 Normal excitation modes 

The harmonic oscillator solutions provide a series of fre- 
quencies after the total separation of coordinates. They 
describe which modes are easiest to excite. For TV molecules 
we always find A'^ — 1 frequencies related to relative mo- 
tion (some of which can be degenerate), and one degree 
of freedom corresponding to the trivial free center of mass 
motion. Note that we consider the limit where the layer is 
in the strict 2D limit, i.e. zero width. This means that all 
modes are transverse. In a real experiment with a finite 
width, longitudinal modes are also possible but usually 
have higher excitation frequencies (in the ID case this has 
been discussed in Ref. [72]). 
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Fig. 11. Frequencies of the six- molecule system as a function 
of dipole strength. The different lines connect frequencies of 
the same normal mode at different dipole strengths. 



In figure [9j we demonstrate how the normal modes 
change with chain length for a given interaction strength. 
As the number of molecules is increased, additional fre- 
quencies appear. When adding one molecule to a system 
of N molecules, the A^— 1 frequencies from this A^-molecule 
system all decrease in magnitude. The new additional fre- 
quency comes in above all the other frequencies, in such a 
way that the average frequency decreases overall. The av- 
erage frequency multiplied by A^ — 1 is in fact the ground 
state energy of the A^-body system, except for the correc- 
tion from the shift in equation (H, which to be precise 
should be multiplied by the number of interacting pairs. 
This is indicated by the horizontal line in figure |9] 

The normal modes in Fig [9] are similar to the acoustic 
branch modes in a crystal with a single atom in the ba- 
sis. In typical textbook presentations of acoustic phonons, 
nearest neighbour interactions are taken into account and 
the spectrum has the form u;(/c)^ = ^sin(^), where 
K is the harmonic force constant of Hooke's law. Taking 
the naive harmonic approximation to the dipolar poten- 
tial (obtained by expanding equation ([I]) to second order 
in r/d), we obtain K = 120^/(1^. The speed of sound then 

becomes cq = \f^^§T as discussed in Ref. (55]. By fitting 
the lowest modes of the chain as function of the chain 
length, we can get an estimate of the speed of sound in 
the harmonic framework used here. 

The possible wave vector in a string of length Nd is A; = 
mr/Nd, where n is an integer n > 0. We can now use the 
fit in figure 10 to obtain the speed of sound and we find c = 
Co/3.05, a reduction by a factor of around three. The first 
obvious reason for the difference is that we have taken not 
only nearest neighbour interactions but all interactions 
into account in our model. However, as we saw in relation 
to the energetics, the terms beyond nearest neighbour play 
a minor role which was also noted in Ref. |62j . The main 
difference is that the effective oscillator potential used here 
allows molecules to spread further in space, it is in a sense 
more shallow than the naive approximation around r/d ^ 
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Mode 

Fig. 12. Normal modes for a chain with six molecules and U — 
10. The vertical axis shows the molecule's layer position. In 
the horizontal direction, the modes are arranged in increasing 
order of frequency (not to scale). The dashed vertical line shows 
the equilibrium position of the molecules, which are shown as 
points. The frequencies of the modes are, from left to right, 
2.00, 3.77, 5.23, 6.34, and 7.03 in units oih/(m(f). 



1. This results in a less rigid string with smaller tension 
and therefore a reduction in the speed of sound. However, 
from the energetics alone it seems clear that this is a much 
better approximation, and we have thus shown that also 
the single-chain modes should be re-evaluated. Notice also 
that in the naive approach, a;(fc)^ scales with U , whereas 
we find a scaling that is softer (C/" with < a < 1) due 
to the reduced oscillator frequency. 

The tension of the chains/strings is directly related to 
the average radius (r^) of the particular potential used 
since this indicates the strength of the interaction and 
how fast the chain returns to its equilibrium position after 
a disturbance. As shown in Ref. [30 , the real potential 
has (r^) slightly smaller than the harmonic approximation 
used here. The naive harmonic approach has (r^) that is 
an order of magnitude smaller than both the real potential 
and the effective harmonic potential used here. This is 
another indication that the spectrum produced here is a 
significant improvement over the naive approach. Based 
on the fact that the real dipole potential has slightly lower 
(r^), the true speed of sound should also be somewhat 
higher than our estimate. 

In figurc[Tl]we show how the dipole strength affects the 
normal mode frequencies. They all increase, as an increase 
of the dipole strength essentially makes the chain stiffer. 
On the other hand, the curves in figure [TT] are not parallel, 
and the frequencies diverge away from each other. The 
higher the frequency of the mode, the more it grows with 
increasing dipole strength. This is different from the usual 
model of acoustic phonons with only nearest neighbours, 
and is due to the inclusion of next-nearest neighbours and 
so on. The latter interactions tend to suppress short wave 
length modes. 

The structure of the normal mode excitations are dis- 
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Fig. 13. Evolution of the normal mode frequencies as a func- 
tion of the intra-layer repulsive frequency, uir- The modes are 
for two strings of six molecules with (7 = 20. 



the maximum displacement amplitudes in one direction 
of the 2D structures. The lowest energy mode has one 
node in the vibration of the string, the next mode has two 
nodes, and so forth. This ordering remains for all dipole 
strengths, and the amplitudes of the molecules change 
very little. 

It should be noted that this pictorial representation is 
in one dimension. The other direction is degenerate but in- 
dependent, which implies that linear combinations of each 
of these one-dimensional modes are equally appropriate. 
The displacements could happen along any two indepen- 
dent directions in the xy-plane. This is reflected in the fact 
that we always find two degenerate modes for each nor- 
mal mode frequency, both for single chains and for several 
chains per layer as discussed below. For example, we could 
choose a radial and tangential mode in the planes instead 
of the Cartesian one-dimensional modes illustrated in fig- 
ure 



12 Still they would have to present nodes at the same 



layers, where the molecules remain at their equilibrium 
position. 



4 Multiple chains and layers 

Molecules in different layers still attract each other as 
parametrized in equation ([3|. However, more than one 
molecule in the same layer immediately introduces the 
repulsion between these molecules. We parametrize this 
additional interaction in equation ([6]) where we treat the 
repulsive frequency, ujr^ as a free parameter. All structure 
calculations are independent of the constant shift which 
only enters in energy calculations. We therefore first dis- 
cuss structure variation as uJr is increased from zero. Even- 
tually there is a point when the repulsion is too strong 
compared to the inter-layer attraction, and the system is 
no longer bound. This is seen as at least one imaginary 
normal mode frequency solution to the full A^-body prob- 



played in figure 12 for a chain of six molecules. We show lem. 
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4.1 Energies and repulsion 

We now solve for N layers and two dipolar molecules in 
each layer. The resulting frequencies are shown in figure 
[13] for A'^ = 6 as function of the repulsive frequency. There 
are 11 normal mode frequencies for these 12 molecules, 
one mode is related to the free center of mass motion. 
Five of these modes do not depend on the repulsive fre- 
quency. They correspond to the modes of the singly occu- 
pied layers (similar to those seen in figure 12), but with 
the frequencies increased by a factor of •\/2rThis reflects 
that the single string degrees of freedom is repeated by the 
pairs in each layer. The strings move in phase as if they 
were alone, and are like the modes shown in figure [T2] with 
pairs following each other instead of single molecules. The 
factor, v^, in the energy arises since the total energy is 
the square root of the sum of energies [57] . Equivalently, 
for two chains the number of attractive pairs double and 
in turn the string tension doubles. The factor comes from 
the acoustic mode frequency which depends on the square 
root of the tension. For M chains, the enhancement is 
then \/M. While this seems at odds with a classical pic- 
ture of a collection of chains in 2D layers, the molecules 
and chains in our model are not localized and therefore 
adding a chain enhances the chain in the same way, irre- 
spective of the number of chains one starts from. 

The remaining normal modes in figure [13] appear as 
three doubly degenerate modes decreasing with the repul- 
sion. They represent relative motion where the two single 
string degrees of freedom are mixed. The degeneracy is 
due to the reflection symmetry in a central plane parallel 
to the layers. Two of these three degenerate frequencies 
coincide resulting in quadruple degeneracy. These modes 
correspond to motion of the molecules in the four middle 
layers. The last degenerate frequency in figure[T3]decreases 
to become the lowest mode, and eventually reach zero for 
a critical repulsive frequency Wc- The system afterwards 
is unstable corresponding to imaginary solutions of the 
energy. 

The behaviour of the frequencies with increasing the 
number of layers is seen in figure 14 (similar to figure [9]) . 
Unlike the single string cases, we see that they are almost 
completely flat. This is analagous to the optical frequency 
modes in crystals with more than one atom at each lat- 
tice site, and we may refer to these as optical (normal) 
modes. As the number of layers is increased, more and 
more frequencies appear in the upper branch, which are 
not resolved on this figure's vertical scale. They are nearly 
degenerate in any case since they are modes of the central 
layers. The effect of repulsion is to push the frequencies of 
the outermost layers to zero, which was also seen in figure 

m 

The structures of the two optical modes correspond- 
ing to the degenerate energy that approaches zero are rel- 
ative motions between the two molecules in each of the 
outermost layers. The instability is then that the ampli- 
tudes of the vibrational motion between the pairs of outer 
molecules become too large to return to equilibrium. The 
degeneracy means that all four molecules then separate 



12 



IS 

■D 























X 


X 


X 


X 

C0r/(i)j,=.99 


X 


+ 


+ 


+ 


+ 


+ 



10 



15 20 
N 



25 



30 



Fig. 14. Normal mode frequencies (points) as a function of 
number of layers for two strings with a dipole strength U = 20 
and two different repulsive strengths. The modes shown here 
are only the intra layer modes. The cohearant full string modes 
behave in the same way as shown in figure [9] 




simultaneously as illustrated in figure 15 Then the whole 



Fig. 15. Schematic illustration of the unstable mode for three 
chains with four particles in each chain. The picture shows 
a ID setup for clarity (the 2D case has a similar cylindrical 
structure), a) The instability of the outermost molecules in 
the top and bottom layers, b) The final state of the system 
after dissociating the two outer molcules. 



system immediately falls apart, as all shorter chains break 
at smaller repulsive frequencies. 

If we increase the number of layers the results would 
be very similar. Only two differences appear. The first 
is that the number of repulsion-independent (horizontal) 
solutions increase, because the number of single chain so- 
lutions increase. The second difference is that the degener- 
acy increases for the highest of the optical normal modes, 
which involve molecules in the intermediate layers. They 
are essentially independent of the length of the chains, 
see figure |14[ The lowest frequency optical mode which 
eventually approaches zero is also almost independent of 
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chain length, and the unstable mode is again related to 
the outermost molecules. 

If we instead increase the number of chains, the fre- 
quencies of the resulting solutions would appear almost 
precisely as in figure [13} Again only two differences appear. 
The first is that now the repulsive-independent (horizon- 
tal) solutions are obtained from single chain frequencies by 
multiplying by the square root of the number of chains. 
The reason is again that energies are obtained as square 
roots of average frequencies as explained above. The sec- 
ond difference is that the degeneracies increase for the two 
decreasing optical mode solutions. The highest of these 
correspond to normal modes involving molecules in the in- 
termediate layers. They are essentially independent of the 
number of chains since they move pairwise without influ- 
ence from the other layers. For each additional chain the 
degeneracy of the terminal layers increases by two, that is 
one corresponding to each end-layer structure. The unsta- 
ble modes are now more complicated and at least one of 
them always involves at least one molecule of small am- 
plitude. This implies that at least one molecule remains 
after the amplification of the amplitudes of the unstable 
mode and a subsequent breakup of the complex. 



4.2 Critical stability 

The structure close to the instability can be seen directly 
from the behaviour of ((r^— i?)^) for the individual molecules 
as function of the repulsive frequency. Here R is the cen- 
ter of mass coordinate projected on a plane. This quan- 
tity can be obtained since each molecule is treated as 
distinguishable. Still, identical bosons in the same layer 
always have the same radii. As a;^ increases all radii of 
interior molecules increase marginally whereas the termi- 
nal positions increase substantially faster. However, when 
uJr is very close to the critical frequency, lOc, the internal 
molecule-radii perhaps double their size but the terminal 
sizes meanwhile increase by an order of magnitude and 
eventually diverge. The reflection symmetry in a central 
plane can be maintained throughout. Due to the circu- 
lar symmetry, the system can be described as a cylinder, 
whose radius changes down its length, becoming wider to- 
wards the ends. 

The critical frequency, ujc, can be found numerically 
by solving the Schrodinger equation. For N layers and 
the same number of molecules in each layer, we label the 
molecules from 1 to iV along the first chain and continue 
the consecutive labeling chain after chain. A formula for 
the critical frequency, uJc, and the shift, a, can be obtained 
analytically by equating the repulsive and attractive po- 
tential energies as seen by one of the molecules in the 
outermost layer where the instability takes place 



where ujij are the input frequencies reflecting the interac- 
tion between a particle in the outermost layer and layer 
j. Eliminating common factors allows one to obtain an 
expression for the critical frequency and the shift factor 



N 

E 

3=2 



and 
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(16) 



(17) 
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where the frequencies on the right hand side of equation 
( |16[ ) are the input frequencies determined in equation Q . 
These expressions hold for a system of equal masses and 
would have to be changed for unequal masses. The result 
provides the value for lUc{N) which explicitly is labeled 
with the dependence on layers. This is a remarkably sim- 
ple and very general result independent of the number of 
chains. The expression for a depends only very weakly on 
U and M, and we find ~ 2.26 numerically for most 
values of TV, M, and U . The expression for lOc exactly re- 
produces the numerically determined values. This is com- 
pletely independent of a since, as we will discuss below, 
a is determined by equating energies of different configu- 
rations, i.e. it is a pure shift of the energy which will not 
influence the stability of the system. Whereas for Uc this 
is just a good initial guess at the proper value, which will 
be commented on in following sections. 

Thus, LJc{N) increases with number of layers, but only 
slowly, since smaller and smaller frequencies from more 
distant molecules are added to the sum in equation (16). 



The hierarchy of stability is now easily established, be- 
cause the sequence follows the size of udN). If a given 
number of layers becomes unstable, the structure with 
one layer less is also unstable. On the other hand, we find 
that the critical repulsion is independent of the number of 
chains. 

So far we considered systems with the same number 
of molecules in all layers. The normal modes are then re- 
lated to symmetric structures essentially localized as mo- 
tion within planes and relative motion between groups 
of molecules in different planes. A different system could 
have fewer molecules in the outer layers than the series 
of identical layers in the interior, an example is shown 
in figure ITsb). The hierarchy of stability follows from the 
number of repulsions in the outer layers, that means one 
molecule is most stable, two second most, etc. 

The critical frequency for a configuration with one 
molecule in the outer layers, uic, is found to be 



1 



^c(^^-2) + — (C^f2+^2JV) 



(18) 



;^^^3iM 



If^ij^li^lj + y%){M - 1) + \^l^jUJla^d 



where we used the definition in equation ( 16 ), the labeling 



N 



3=2 



N 
3=2 



^ ^ are that molecules 1 and N are the lone outer molecules, 

and molecule 2 is in the second layer with M molecules, 
(j ipg^restingly, the ratio, uj-^/uj^{N), seems to be indepen- 
dent of the dipole strength. 
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Fig. 16. Single chain modes as a function of number of chains, 
M, for a broken system of 6 layers (where now there is only 
one molecule in the outermost layers). The lines connect the 
same normal mode at a given number of chains. The dashed 
black line is the frequency of the degenerate mode which is not 
the breaking mode frequency taken when cjr ~ (jJc- The dipole 
strength was U — 20. 



The behaviour of the normal modes and normal mode 
frequencies is similar to the complete system. The total 
number of remaining modes is M{N — 2) + 1. The — 1 
single chain modes remain, and they appear qualitatively 



just as seen in figure 12 These frequencies increase with 
the number of chains, and are shown in figure |16[ but 
no longer in the square root dependence that was seen in 
the complete system. As with the unbroken system, the 
single chain frequencies do not change with the repulsive 
frequency. The remaining modes fall into two degenerate 
groups of equal size, 2(M— 1). One group approaches zero 
as the repulsive frequency is increased, while the remain- 
ing modes are at a frequency between the two lowest single 
chain modes. The group approaching zero are modes in- 
volving the outermost full layers, which become unstable 
as the frequency is increased. The other group are modes 
in the central layers. 

Another question is which structure the unstable mode 
has when only one molecule occupies the outer layers. The 
unstable mode is moved to the next-to-outermost layers 
where the molecules close to the critical frequency then 
behave as in the system where they were in the outer 
layers. The resulting configuration is then most likely one 
molecule in each of the two outer layers in each end. This 
configuration is stable because the decrease of io^{N — 
2) to the value uj'^{N — 4) is compensated by the added 



dominating term, ujf2, as seen in equation ( 18 ). Increasing 



the repulsion would repeat the process of instability of the 
next fully occupied layers. The system seems to prefer to 
approach single chain structure. 

If we decrease the number of molecules in the outer lay- 
ers by one molecule at a time, the instability is seen as de- 
generate modes of varying frequency approaching zero for 
sufRciently strong repulsion. The degeneracy corresponds 
to an instability in the outer layer until the outer layer 
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M 

Fig. 17. Density of the central molecules as a function of rmm- 
ber of chains. Each chain is eight molecules long and the dipole 
strength was U = 20. 



has just more than half as many molecules as the second 
layer. Then the degeneracy changes to correspond to in- 
stability of the molecules in the second layer. This reflects 
that the instabilty is initially related to the outer layer 
until a number of molecules is removed from that layer 
and the second layer contains the unstable mode. Regard- 
ing the critical frequencies, for the outermost layer, the 
critical frequency is 



h 



M -h 



k=2 



(19) 



where h is the number of holes in the outermost layer, 
and ujij is the interaction frequency between layer one 
and layer j. For the molecules in the penultimate layer, 
the critical frequency is 



M -h 
M 



^2N 



(20) 



Whichever of these two frequencies is lower for a given 
configuration determines in which layer the instability lies. 

In summary, for the frequencies, there are a total of 
NM—1 normal mode frequencies. A subset, A^— 1 of them, 
are single string modes with frequency equal to ^/M the 
frequency of the corresponding single string mode. There 
are then many degenerate modes, 2(M — 1) of them go 
to zero as the repulsive frequency is increased. This is 
because the chain breaks by removing the repulsion in 
the outermost layers, releasing all but one molecules. The 
remaining degenerate frequencies cluster into groups, and 
are often very close to each other in magnitude. They also 
decrease as the repulsive frequency is increased, similar to 
what is seen in figure [T3| for two strings. 



5 Density and energy relations 

The repulsive frequency and the repulsive shift are so far 
left as free parameters, although we did discuss at length 
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how one can relate them to two-body properties in section 
[2] The shift is only related to total energies without any 
influence on structures. We need to define a zero point to 
measure from but this requires more than two molecules, 
and the criterion easily depends on which system was used 
for the calibration. In the next subsections we deal with 
these two parameters in the in-layer repulsion. 



5.1 Density versus repulsive frequency 

In a system with many layers, the density depends on the 
relative position of the layer. To be specific we consider the 
central layers where the density is rather constant across 
a few layers. The outer layers close to the critical stability 
are on the other hand expanding as the breaking point 
is approached. Already for a relatively small number of 
layers this expansion only affects the outer layers. We thus 
define the density through the average radii of the central 
layers, i.e. 



M 



M 



(21) 



fc=r 



where this definition is valid for non-identical molecules. 



The densities of the central layers defined in equation (21 1 



are seen in figure [17] for different repulsions as functions 
of number of molecules in each layer. These curves are 
essentially independent of chain length, but increase with 
molecule number and decreasing repulsion. 

The density of the central layers are essentially inde- 
pendent of chain length. It does depend, however, on the 
dipole strength, molecules per layer, M, and the repulsive 
frequency, i.e., p oc a{M,u}r, U) with a very weak depen- 
dence on TV. The dependence on M is a straightforward 
M3/2, which comes from, in eq j2lL with one power of 



M in the numerator, and it is found that (r^) goes like 
^^-1/2 rj^Yie critical frequency, and hence the critical den- 
sity, below which the system becomes unstable, depends 
on the parameters as mentioned above. We can state ex- 
plicitly the dependence on Af , and thus try to isolate its 
dependence on frequency, and thus U: 



cx Af3/V(iV, C/), 



(22) 



though the dependence on N is very weak. The depen- 
dence on ujc is linear, and the dependence of LodN, U) on 
U is seen in figure 18 One can see that the dependence on 
the length of the chain, N, is indeed very weak. However, 
in looking at densities of central layers, one must have 
a chain of 4-5 layers in order for a "central layer" to be 
clearly defined. A plot of the density as a function of inter- 
mediate values of the repulsion is seen in figure[T9] One can 
see that when the chain lengths are doubled, the curves 
are almost indistinguishable. The relations between the 
lines that differ in the other parameters are as described 
in eq ( 22 ) . The behaviour of the curve is not analytical, 
but the ratio p{LOr = 0)/Pc = 1.39, where p{uJr = 0) is 
the density with zero repulsion, is seen for all the sets of 
parameters. 
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Fig. 18. Critical frequency for two chains of differing lengths 
as a function of interaction strength. 



In section [2j we discussed the option of relating the 
magnitude of the repulsive frequency to density (see ([s])), 
but did, however, caution that this is not a very desirable 
approach to effective interactions in the iV-body problem. 
The density that we obtain within our model can become 
very large. In fact, if we insert the obtained densities of or- 
der crd^ ~ 10— 100 in Q, the frequency obtained is larger 
than the critical frequency by an order of magnitude al- 
ready for M = 3 and increases rapidly with M. This shows 
the dominance of the attractive interlayer forces, i.e. the 
system contracts into denser configuration than the typ- 
ical average experimental density. We stress again that 
([s]) is a questionable relation but note that the alterna- 
tive way of fixing the repulsion through the bound state 
structure of the potentials also predicts a repulsion slightly 
beyond the critical value. So in a system where the inter- 
and intralayer dipolar strength is the same, we expect the 
configurations with multiple molecules in single layer to 
be unstable. 

An important point here is that the critical frequency 
is essentially independent of M, i.e., if we add another 
chain to the system for a fixed value of uJr < uJc{N, U). 
On the other hand, the density goes up since the addi- 
tional attractive interlayer terms that balance the extra 
repulsion in each layer will make the system contract. Ac- 
cording to ([s]) this implies that should be increased and 
it will then quickly exceed u)c{N, U). Again, it is unfortu- 
nate to change the two-body interactions in response to 
a property of the total system. By allowing the repulsion 
to be a free parameter we can study the detailed competi- 
tion between repulsive and attractive term in the system 
which was one goal of the current study. 



5.2 Energy versus repulsive shift 

The relative energies depend on the repulsive shift in equa- 
tion ([7]). A consistent choice of scale factor, a, independent 
of chain M, and layer, N , numbers is not possible. The 
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Fig. 19. Densities as a function of the ratio of the repulsive 
frequency to the critical repulsive frequency. This is plotted for 
several values of the parameters A'', Af, and U, and the curves 
are identified by the label {N, M, U). 



problems arise, because instability towards some contin- 
uum structure is found through the lack of solutions, and 
totally independent of the energy of the configurations. 
However, the repulsive shift is precisely able to adjust en- 
ergies such that the wave function and energy criterion 
coincide for specific geometries. The main effects can be 
included by an appropriate choice. 

The constant, a, may be obtained by comparing the 
energies of one system at the point of instability with the 
other system formed after breaking apart. The constant 
is adjusted so that the energies of these systems are the 
same at the breaking frequency. Let us choose the simplest 
structure with three molecules in each layer and compare 
to the system where two molecules are removed from each 
of the end-layers. This requirement leads to a ~ 1.92, 
which makes the repulsive shift somewhat weak compared 
to the attractive one at the breaking frequency. This choice 
of a has the nice feature that as the number of molecules 
per layer is increased, the energy becomes more negative, 
but does not run away to negative infinity. This is a nice 
feature since we know that the breaking frequency is inde- 
pendent of M . However, this value of a will not match en- 
ergies when the system is driven to its breaking frequency. 
A better form for the constant might be 



= ao ± al/M, 



(23) 



stant do what is prescribed of them, but do have draw- 
backs. For example, for M = 2, the overall sign is negative, 
which makes it an additional attractive shift. Also, while 
the energies at breaking match at M = 8, at A/ = 9 they 
already are mismatched by 1.2%. Clearly one value of a 
cannot accurately reproduce all the desired properties of 
the total energies. A reasonable average choice is a ~ 1.92 
but more accuracy can be achieved by focussing locally on 
a specific system and tune the value of a. 

This value of a ^ 1.92 was cited in the discussion of 
the two-body repulsive interactions in section [2] and leads 
to a reasonably low w^, which is, however, still larger than 
the critical frequency for breakup. Again we remark that a 
is a geometric quantity related to the form of the inverted 
oscillator that provides the repulsion. 

By increasing the number of chains, or increasing the 
number of molecules in each layer of the cylinder, we 
are essentially adding another dimension to the system, a 
number of molecules per layer in addition to the length of 
the chains. For a system of M chains of N molecules, one 
can count the number of attractive and repulsive pairs in 
the system. If we keep to the nearest neighbour level of at- 
traction, then the number of attractive pairs is M'^{N — 1), 
and the number of repulsive terms is M{M— l)N/2. Their 
ratio is then 



Attractive pairs/Repulsive pairs 



2M{N - 1) 
N{M - 1) 



(24) 



where ag is obtained by fitting to the quadratic depen- 
dence of the energies of M in systems where M > 5. 
The second constant, ai, is determined by matching the 
energies of a complete system, and the system where one 
molecule has been removed from both of the outermost 
layers (where it is assumed that the ao term will can- 
cel). This combines two desirable features, controlling the 
energy dependence on M and matching with the broken 
system. These constants were determined for the system 
{N, M) = (6, 8) and [/ = 20 to be (ao, ai) = (1.9761, 3.0123) 
and the lower sign is taken in equation (23). These con- 



This ratio appears, for example, in the bilayer system 
{N = 2) , where if one keeps adding more pairs of molecules 
(one in each layer per pair) , the ratio of attraction to re- 
pulsion decreases from two (M = 2) to one (M ^ 1). One 
might expect that this would influence the critical repul- 
sive frequency, but this frequency does not change with the 
addition of more chains. In the nearest neighbour approxi- 
mation, it would also not change with increasing the chain 
length. It is the interaction beyond the nearest neighbour 
that can affect the critical frequency. The relation can, 
however, affect the energetics of the system, depending 
on the relative strengths of the attraction and repulsion. 
As mentioned before, this indicates that a scale factor de- 
fined for a system with only three chains (M = 3), would 
be too large if used for systems with large M, as the ratio 
of attractions to repulsions falls as M increases (for fixed 
N). 

An example application of the scale factor is for a sys- 
tem of ten chains of six molecules, where molecules are 
removed as the repulsion is increased. The molecules are 
removed from the layer indicated by the degeneracy of 
the normal mode frequency that is going to zero at the 
critical repulsive frequency. It was found that the energet- 
ics correctly describe the order of the removal, which can 
be seen in figure [20j When plotted relative to the ground 
state energy (and scaled down because the energies are 
quite large in our units), the point where the degenera- 
cies switched to the penultimate layer was also energeti- 
cally lower than continuing to remove molecules from the 
,outermost layer. For the figure, the constants in the re- 
pulsive shift were those quoted above for the system of 
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Fig. 20. Energies relative to the full six layers often molecule 
system, i.e., {N,M) = (6,10) as molecules are removed from 
the outermost layers, molecules are removed one-by-one from 
both outer layers, in order to preserve the symmetry about the 
center layer of the system. When performing the calculation in 
a system where molecules have been removed, the repulsive 
frequency used is the critical frequency of the parent system. 
This was not done for the the system where ten molecules 
where removed, since the uJcifi holes) > udlQ holes). For this 
system with ten holes, the calculations were done for the two 
different geometries: five holes in each outer layer, and with 
four holes in the outer layer and one in the penultimate layer, 
at the critical frequency of the latter geometry, which is the 
lower of the two and lower than that of the parent system. 

{N,M) — (6,8), so one can see the mismatch for the re- 
moved molecules aheady appearing, as there should be no 
energy difference between the full system and the system 
where one molecule has been removed from the outermost 
layer at the critical frequency. Thus if one is interested in 
detailed energetics, one should decide on a shift energy in 
the region where one is interested in doing calculations. 



6 Summary and Conclusion 

We have discussed how the harmonic approximation can 
be used to study the iV-body 2D problem of dipole-dipole 
interacting molecules. First the two-body interaction is 
approximated by a shifted harmonic oscillator potential 
with the same negative and positive spatial regions. The 
frequency is then determined as function of the molecular 
dipole moment by the requirement that the two-body en- 
ergy has to be reproduced. We test this procedure by com- 
paring with the exact energy of three and four molecules 
in separate layers. The approximation is remarkably ac- 
curate over the range of available dipole moments. The 
repulsion between in-layer molecules is also approximated 
by a shifted harmonic oscillator but now we leave the fre- 
quency as a parameter which can be related to properties 
of larger bound state structures when they become avail- 
able from exact calculations. 



The basic ingredient is a chain where one molecule is 
placed in each layer where the intralayer repulsion then 
is absent. We calculate the binding energy in units of the 
two-body energy of the number of neighbouring pairs, and 
find a weak increase as function of molecule number. The 
value is always only slightly larger than the analytical es- 
timate of C(3) which is approached for large dipole mo- 
ment. Weaker interactions result in about 10% larger val- 
ues reflecting more correlation along the chain. The near- 
est neighbours contribute by more than 80% of the energy. 
The spatial extension in the planes measured in units of 
the layer distance from the total center of mass position 
increases from central towards the outer layer. The weaker 
the interaction the larger is the size. The normal modes are 
oscillations away from equilibrium with a different num- 
ber of nodes. All frequencies of the normal modes increase 
with interaction strength, while the lowest of these always 
decrease with increasing molecule number. 

The normal modes of the single chain configuration 
are equivalent to those obtained in the theory of acous- 
tic phonons in crystals. By fitting our numerical results 
for the lowest mode frequencies, one can obtain the speed 
of sound in the chain. We found that the speed of sound 
was more than a factor of three less than what one would 
naively get by expanding the dipole potential around the 
origin to quadratic order. This reflects one of the dif- 
ferences of our procedure for obtaining an effective har- 
monic approximation through the exact two-body bound 
state energy and the potential proflle. The latter approach 
gives an effective more shallow potential which results in a 
smaller string tension and a reduced speed of sound. This 
indicates a short-coming of the naive harmonic approach 
in describing the collective behaviour of a chain of dipoles. 

We proceed to investigate two chains consisting of two 
molecules in each layer. The structure of the system is 
most clearly found through its excitation modes, or specif- 
ically the normal modes. The repulsive intralayer interac- 
tion is varied from zero to the critical value where the sys- 
tem becomes unstable as seen by a vanishing normal mode 
frequency. The normal modes fall into two distinct groups. 
The first group is independent of the intralayer repulsion 
with a simple one-to-one correspondence to the one-chain 
modes. These modes correspond to in-phase motion of the 
two chains. The energies of the second group decrease with 
repulsion and the lowest of these modes vanishes at the 
critical value. These modes correspond to in-layer oscilla- 
tions breaking the one-chain structure. This is analagous 
to optical modes in crystals with more than one atoms at 
each lattice site. The breakup is seen through diverging 
amplitudes on the molecules in the two outer layers. 

Increasing the number of chains so that there are more 
molecules in each layer leads to remarkably similar be- 
haviour of the normal mode frequencies. The one-to-one 
correspondence with one-chain modes remains. Further- 
more, the decreasing normal mode frequencies also van- 
ish at the same critical repulsion which then is indepen- 
dent of the the number of molecules in the layers. How- 
ever, the degeneracy increases in correspondence to the 
increase in the number of degrees of freedom within one 
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layer. This implies that the configuration at the break- 
ing point has diverging amplitudes of the molecules in 
the two outer layers. The breaking then leads to a sys- 
tem where the two outer layers have fewer molecules than 
the more central layers. These configurations have less re- 
pulsion and the next instability is related to removal of 
molecules from the next two outer layers. This peeling off 
molecules would continue to the third outermost layer or 
remove more molecules in the outer layers. Eventually one 
chain only would remain. We also calculate the density of 
the central layers and find that it is independent of the 
number of layers, depending on the number of molecules 
per layer to the power 3/2 with a factor that depends on 
the dipolar strength. 

The structures are independent of the constant shifts 
in the repulsive potential terms which only have influence 
on the relative energies. We choose these shifts such that 
the energy of the system at the critical frequency is the 
same as the most likely combination of fragments. We ad- 
just to relatively large numbers of layers and molecules in 
each layer. Then we compare energies of a system and two 
systems of half the sizes either in layers or molecules in the 
layers. We conclude that the more layers, the more stable 
is the system, whereas the energy is rather independent of 
the number of molecules in each layer. 

In conclusion, the harmonic approximation is very use- 
ful in the description of the two-dimensional systems of 
cold polar molecules. However, this requires a careful ad- 
justment of the forces between molecules in different planes 
and in the same plane. The latter can be particularly 
tricky to define and we have discussed at length how to de- 
scribe the in-plane repulsion in a reasonable way using an 
inverted harmonic oscillator potential. Using a harmonic 
approximation for both types of interaction is very desir- 
able due to its exact solvability. The structures are found 
as function of layers and molecules in the layers, energies 
and sizes are calculated, and configurations of excitation 
modes are found. The structures are ordered in a hierarchy 
related to their stability. 

The results obtained here within the harmonic model 
indicate that for a typical system where the strength of 
the inter-layer and intra-layer terms are both proportional 
to the square of the dipole moment, bound complexes 
with more than one molecules in a single layer are not 
bound and that the single chains with one molecule in 
each layer are the stable entities. However, the external 
confinement that any realistic experiment employs means 
that the dipoles cannot separate to large distances. We 
therefore expect that the multi-chain modes calculated in 
the harmonic model could show up as resonances. Here we 
provide a study of these modes and the breakup channels 
of multi-chain complexes. 

The experimental conditions needed to study the sys- 
tems considered here should be within reach of the next 
generation of cold polar molecular experiments. The cur- 
rent front-runner uses '^'^K^^'Rb molecules which have U = 
0.1 for the condition in Ref. lOj, but which can probably 
be extended up to about U = 1.2 with increasing electric 
field. These values of U are likely too small for the har- 



monic approximation used here to apply since the two- 
body bound states are very extended [27j . However, most 
other combinations of alkali atoms into polar molecules 
can potentially be used at dipole moments of 1 Debye or 
more |73| . which would easily increase U by one or two 
orders of magnitude. Larger values of U implies that the 
single chain structures will have large binding energies, 
and in turn be stable even in thermal samples. This means 
that one does not necessarily need to be in the nano-Kelvin 
temperature regime to study the chains discussed in the 
current work. 

The current work indicates that the structure of a 
multi-layer system with several particles in each layer is 
most likely that of single chains, i.e. the dipolar chains 
liquid originally proposed in Ref. [T? , but that more com- 
plex states are possible if the intralayer repulsion could be 
somehow tuned independently of the interlayer attraction. 
This could perhaps be obtained by using a combination 
of external electric fields and applied lasers [68..69..70, . 
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